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Abstract 

An adaptive, kink-based path integral formalism is used to calculate the ground state energies of 
the atoms He-Ne. The method uses an adaptive scheme to virtually eliminate the sign difficulties. 
This is done by using a Monte Carlo scheme to identify states that contribute significantly to the 
canonical partition function and then include them in the wavefunctions to calculate the canonical 
averages. The calculations use the 6-31G basis set and obtain both precision and accuracy. 
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I. INTRODUCTION 



The path integral formulation of quantum mechanics offers a variety of advantages for 
studying the electronic and geometric structures of multi-electron systems 0. Chief among 
these are inclusion of finite temperatures (particularly as they affect geometric degrees of 
freedom) and exact inclusion of electron-electron correlation. The application of this method 
to electronic systems has been hindered by the so-called "sign" problem, which results 
from sign of the fermion density matrix, which can be positive or negative and leads to 
large uncertainties in quantities evaluated using statistical methods such as Monte Carlo 
simulations^, §, |, § §, [7], §, g [OJ [IT], Q [JJ, [JJ, 0, [TJ, 0, © 0. We have re- 



cently introduced a "kink-based" path integral approach ||19||, which was demonstrated to 
overcome the sign problem in the 2-D Hubbard model. This approach is complimentary 
to the shifted-contour auxiliary-field Monte Carlo method^, [16], [17|, pLSj] , which uses the 
Hubbard-Stratonovich transformation to combat the sign problem. In this work, we use 
the kink-based formalism to study atomic systems, the next step in studying systems with 
geometric degrees of freedom (such as atomic clusters). 

II. KINK-BASED APPROACH 

In this section, a brief review of the kink-based approach [ 19] is given, with additional 



attention given to the different spin states that are encountered in electronic systems. The 
partition function is written: 

Q = Tr{exp(-/?#)} 

= < a, a\ exp(— (3H)\ a,a > 

CJ.XX 

= J2 e M-PE*,°) (i) 

<7,a 

where a labels the different electronic states associated with a particular spin state a and 
| a, a > is the properly anti-symmetrized state. For large enough (3, this becomes 

Q w exp(-/3£ 0>CT *) (2) 
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where E 0jCr * is the ground state energy of the lowest energy spin-state. If an approximate 
set of states, {a, s} is used, we have 

Q{a,s} = ^2 < °' S ' eX P(~^)l °> S > 

a,s 

= < a,s\a,a > \ 2 exp(-(3E atCT ) (3) 

a,s tr,a 

As long as < a, s|0, a* > ^ for some a and s, then as (5 gets large, 

Q {ajS} oc exp(-/3S 0i(T *) (4) 

In a later section, we will choose our states with specific values of S z . Consequently, we 
will determine the low temperature partition function corresponding to the lowest energy 
spin-state S that has S z as one of its possible values of S z . 

To evaluate the partition function Q{ a , s } using the path integral method, we insert com- 
plete sets of states in order to use the high temperature, semi-classical approximation for 
the density matrix: 

Q{a,s} = <a 1 ,s 1 \ex.p(-pH/P)\a 2 ,s 2 > ■■■ 

ai,si ap,s P 

x < ap, sp\ exp(— @H/P)\ai, si > 

= . . . f a 2, s 2 . . . /m,si (c\ 

a\,s\ ap,sp 

We refer to a matrix element < a, s\ exp(—(3H/P)\a', s' > with a ^ a' or s ^ s' as a kink. 
We rewrite the partition function as a sum over kinks: 

a,s 

EEEte) , («)" 2+ *( t «) 2 + '" < 6 > 

j=l a,s a' ,s' 

= Qo + Q2 + Qs + ■ ■ ■ + Qp (7) 

where Q n is the partition function corresponding to n kinks. In our previous work, we 
demonstrated that Q{ a ,s} has the form 

Q{a,s} = J2 X 3 + 
j 

(nz) fn s({ Xj },n, m ,{ gj }) (8) 

„=2 \i=l ji J \k=l J 
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where 



Xj =< a>j, Sj\ exp(—j3H/P)\ atj, Sj >~< aj, Sj\(l — (3H/ P)\ a>j, Sj > 
tjji =< aj, Sj\exp(—j3H/P)\ ctj>,Sj> >~< aj,Sj\(l — j3H/P)\ ctj>,Sjf > 

and S {{xj} , n, m, {gj}) is the contribution to the partition function with n kinks, comprised 
of m states atj, each occurring gj times • gj = P — n). The explicit form for S is 



S({xj},n,m,{gj}) = 



1 d 91 



-i 



.'I i 



j^( gi -l)\dxf- 1 Y[ k ^{xi-x k \ 
The derivatives can be evaluated recursively. If we define 



Hi- 



dP 



x 



p-i 



dx l Uk^l ( X l - X k 



\9k 



we can show 



m 

* = E ' 



n-l 

E 



m=0 



n — 1 
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(m) 



-l) m m! 



P- 1 



.r 



m+l 



E 
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A similar manipulation leads to an expression for the energy estimator: 
E est = + 



1=1 



1 / •) 

1 ^ ( gi - 2 



5^ (^- 1)! ^ 

Z=0 vi " ; 7=0 



j 



(m) 



(m) 



(m) 



dp 1 



(9) 



(10) 



(11) 
(12) 

(13) 
(14) 



(15) 

(16) 
(17) 



The expression shown in Eqn. |^ gives the exact value of Q (including electron-electron 
correlation) within the approximations inherent in using a finite basis set and a finite level 
of discretization. The so-called "sign problem" can occur in this and any other discretized 
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version of the path integral problem because any of the matrix elements can be nega- 
tive, resulting in a large variances when evaluating the partition function using simulation 
methods. Our approach to minimizing or eliminating the sign problem has been to apply 
Eqn. [| in an adaptive manner. We first realize that the zero kink contribution to Q has no 
sign problems, since the system is in a single state. With a properly chosen states, Q can be 
obtained with just a few kinks; this significantly reduces the sign problem by reducing the 
statistical error from greater than 100% to a precision adequate for chemical applications. 
A good choice for the states is obtained using a Monte Carlo simulation, in which the dif- 
ferent N-electron states that appear during the simulation are used to update the estimates 
of the ground and excited states. We call this approach an adaptive approach, since the 
Monte Carlo algorithm allows the estimates for the ground and excited state wavefunctions 
to evolve according to the statistical sampling of the different N-electron states. 

We implemented the adaptive scheme in the following way (other methods of are possible). 
An initial set of basis functions (the 6-31G basis set in our calculations) was orthonormal- 
ized and used to create a set of one-electron orbitals. The one-electron Hamiltonian was 
then diagonalized in this basis. Each electron was assigned a spin and the one-electron basis 
functions were combined to form a set of Slater determinants (as described earlier, the low- 
est energy spin state will be projected out by the path integral procedure) that were then 
used as the initial \a, s > for the Monte Carlo simulation. A simulation using the absolute 
value of the summand in Eqn. ^ as the weighting function was performed in which kinks 
were added, removed, and changed. An upper limit on the number of kinks allowed was set 
to 10 kinks; since the final results were obtained with or 2 kinks, this did not affect the 
accuracy of our results. A list of the states accepted was kept. If the fraction of configura- 
tions that contained more than kinks was greater than the fraction of configurations that 
had kinks, the Hamiltonian was diagonalized using the current list of states and a new 
set of N-electron states obtained. These new diagonalized states were linear combinations 
of the initial set of Slater orbitals and thus corresponded to configuration interaction (CI) 
wavefunctions. Since the simulation sums over all possible states, in essence a complete CI 
calculation is performed. Another possible Monte Carlo scheme would allow the individual 
Slater determinants to be altered during the simulation, which would correspond to a MC- 
SCF calculation. At most 100 states were included in the diagonalization to limit the time 
per diagonalization. If the set of accepted states exceeded 100 at the time of diagonalization, 
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only the 100 most prevalent states were included. Once 5000 iterations had occurred with 
no diagonalizations, the run was terminated and these final 5000 energies used to determine 
energies. At the end of the calculation, the ground state would correspond to a high quality 
CI ground state; if the state does not correspond to the complete CI wavefunction, then 
kinks will be added to correct the ground state. If the adaptive procedure provides the com- 
plete CI wavefunction, then no kinks would ever be introduced, as the density matrix would 
be diagonal. In practical calculations, we expect to stop the adaptive process before the 
density matrix is actually diagonal, but when the off-diagonal matrix elements are so small 
that the likelihood of adding more than 2 kinks is very small. In fact, in our calculations, 
the Monte Carlo procedure provided such a good estimate of the true ground state that at 
any time we found only kinks or 2 kinks (to one of the excited states). 



III. APPLICATION TO ATOMIC ENERGIES 



We have tested this approach by applying it to atomic systems, using the 6-31G basis 
set. This set was chosen for its relative simplicity and reasonable accuracy. For each atom, 
He-Ne, each electron was assigned a specific s z , leading to a fixed total S z . Thus, the sum 
in Eqn. |5| used just a single spin state which was a linear combination of states S such that 
S z was one of the possible values of S z . The initial basis functions were orthonormalized 
and the one-particle Hamiltonian was diagonalized, providing an initial set of states. As 
the Monte Carlo simulation was performed and diagonalizations proceeded as previously 
described. Our results are shown in Table |, along with the Hartree-Fock and CASSCF 
energies from Gaussian 98|2(J. The average sign of the density matrix demonstrates that 
the adaptive approach adequately reduces sign problem to well below what is needed for 
chemical accuracy. For comparison, we note that a shifted-contour auxiliary-field Monte 



Carlo calculation of NejfRJ, using a 4-31G basis set, led to errors of 0.004 a.u., significantly 
larger than those found in the present calculations. For illustrative purposes, Table [TT] shows 
the evolution of the coefficients of the Slater determinants that contribute significantly to 
the ground state wavefunction of Be, during the first 4 updates. After the first 4 updates, 
only minor changes occurred in the ground state. It can be seen that the adaptive procedure 
introduces mixing between the Slater determinants as needed. The degeneracies seen can 
be rationalized on the basis of symmetry. 
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Atom 


E(HF) 


E(CASSCF) 


E(MC, P=10 13 ) 


< Sign > 


N up 


■^down 


He 


-2.855160 


-2.8701621 


-2.8701621(0) 


1.0000(0) 


1 


l 


Li 


-7.4312350 


-7.4315542 


-7.4315535(6) 


1.0000(0) 


2 


l 


Be 


-14.5667641 


-14.6135453 


-14.6135468(22) 


1.0000(0) 


2 


2 


B 


-24.5193448 


-24.5628917 


-24.5628918(14) 


1.0000(0) 


3 


2 


C 


-37.6768656 


-37.7162644 


-37.7162663(24) 


1.0000(0) 


4 


2 


N 


-54.3820508 


-54.4199396 


-54.4199404(32) 


1.0000(0) 


5 


2 





-74.7782342 


-74.8394081 


-74.8394091(38) 


0.9992(11) 


5 


3 


F 


-99.3602182 


-99.4474231 


-99.4474225(34) 


0.9996(8) 


5 


4 


Ne 


-128.4738769 


-128.5898023 


-128.5898026(24) 


0.9996(8) 


5 


5 



TABLE I: Hartree-Fock (HF), CASSCF, and path integral (MC) energies (in atomic units), and 
average sign of the density matrix for the different atoms studied in this work. The numbers in 
parenthesis represent 2 standard deviations. The number of up- and down-spin electrons is also 
specified. 



Update Number 


E (ground state) 


State 1 


State 2 


State 3 


State 4 


State 5 




Degeneracy 


1 


2 


1 


3 


6 


1 


-14.3543 


1.0 


0.0 


0.0 


0.0 


0.0 


2 


-14.6053 


0.7788 


-0.3729 


0.1408 


-0.1340 


-.07652 


3 


-14.6131 


0.7782 


-0.3726 


0.1407 


-0.1339 


-.07645 


4 


-14.6135 


0.7715 


-0.3767 


0.1452 


-0.1343 


-.07872 



TABLE II: Energies (in atomic units) and coefficients of the ground state wavefunction, as a 
function of adaptive update, for the first 4 adaptive updates of Be. States are arbitrarily labeled 
and correspond to multiple states with degeneracies as indicated. The Monte Carlo procedure 
identified all degenerate states and mixed them with identical coefficients. 

IV. CONCLUSIONS 

The adaptive, kink-based approach to path integral calculations has been applied to 
atomic systems. As was the case in our previous work, the use of the adaptive approach 
reduced the sign problem to a tolerable level. While we have used an adaptive diagonalization 
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procedure to improve our estimates for the electronic states, this is not an essential ingredient 
in the adaptive approach. For instance, unitary transformations can be sampled as part of 
the Monte Carlo process. In addition, we have not made any simplifying assumptions that 
will be required when treating systems with large numbers of basis functions, such as limiting 
the type of determinant that can contribute to the ground state wavefunction. We note that 
the number of electrons and basis functions used in this study is on the order of that needed 
to study moderately large metal clusters. For example, Na 2 o would require roughly 10 shell- 
orbitals (orbitals centered at the origin of the cluster, in accord with the shell model of the 
electronic structure) and 20 electrons, if pseudopotentials are used for the core electrons. 
Thus, the current method may be applicable to moderately large systems. 

V. ACKNOWLEDGEMENTS 

The computers used in these calculations were purchased with funds from NSF 9977124. 



[1] R. Feynman, A. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill, 1965. 

[2] C. Mak, R. Egger, H. Weber-Gottschick, Phys Rev Lett 81 (1998) 4533-4536. 

[3] C. Mak, R. Egger, J Chem Phys 110 (1999) 12-14. 

[4] R. Egger, L. Muhlbacher, C. Mak, Phys Rev E 61 (2000) 5961-5966. 

[5] S. Miura, S. Okazaki, J Chem Phys 112 (2000) 10116-10124. 

[6] B. Militzer, W. Magro, D. Ceperley, Contrib Plasma Phys 39 (1999) 151-154. 

[7] W. Newman, A. Kuki, J Chem Phys 96 (1992) 1409-1417. 

[8] R. Hall, J Chem Phys 94 (1991) 1312-1316. 

[9] R. Hall, M. Prince, J Chem Phys 95 (1991) 5999-6004. 

[10] D. Ceperley, Phys Rev Lett 69 (1992) 331-334. 

[11] P. Roy, S. Jang, G. Voth, J Chem Phys 111 (1999) 5303-5305. 

[12] P. Roy, G. Voth, J Chem Phys 110 (1999) 3647-3652. 

[13] B. Militzer, W. Magro, D. Ceperley, Contrib Plasma Phys 89 (1999) 151-154. 

[14] D. Ceperley, in: K. Binder, G. Ciccotti (Eds.), Monte Carlo and molecular dynamics of 
condensed matter systems, 1996. 



8 



[15] N. Rom, E. Fattal, A. K. Gupta, E. A. Carter, D. Neuhauser, J. Chem. Phys. 109 (1998) 
8241-8248. 

[16] N. Rom, D. M. Charutz, D. Neuhauser, Chem. Phys. Lett. 270 (1997) 382-386. 

[17] Y. Asai, Phys. Rev. B 62 (2000) 10674-10679. 

[18] R. Baer, M. Head-Gordon, D. Neuhauser, J. Chem. Phys. 109 (1998) 6219-6226. 

[19] R. Hall, J. Chem. Phys. 116 (2002) 1-7. 

[20] M. J. Frisch, et al., Gaussian 98 (revision a.x) (1998). 



9 



